function output = rhythm_cpu(y);
if nargin < 1, y = wavread('technobeat.wav'); end

bandlimits=[0 200 400 800 1600 3200];
numBands = length(bandlimits);
maxsigFreq=44100;

hannWinLength = .4;


hannlen = hannWinLength*2*maxsigFreq;

yFFTed = fft(y);
n = length(yFFTed);
bandRanges = zeros(12,1);

i = [];

for i = 1:numBands-1
    bandRanges(2*(i)-1) = floor(bandlimits(i)/maxsigFreq*n/2)+1;
    bandRanges(2*i) = floor(bandlimits(i+1)/maxsigFreq*n/2);
    
end

% the end cases of the band limits are special cases
bandRanges(2*numBands-1) = floor(bandlimits(numBands)/maxsigFreq*n/2)+1;
bandRanges(numBands*2) = floor(n/2);

output = zeros(n,numBands); %for storing the output
for i = 1:numBands
    start = bandRanges(2*i-1);
    stop = bandRanges(2*i);
    output(start:stop,i) = yFFTed(start:stop);
    output(n+1-stop:n+1-start,i) = yFFTed(n+1-stop:n+1-start);
end

output(1,1)=0;


sigTime = ones(size(output)); % signal in time domain
sigFreq = ones(size(sigTime)); % signal in frequency domain
hann = zeros(n,1); % for the hann 

%for the (relatively) small value of hannlen, it's slower to use the GPU
for a = 1:hannlen
    hann(a) = (cos(a*pi/hannlen/2)).^2;
end


% Back to frequency domain with absolute values in the time domain;
for k = 1:numBands
    sigFreq(:,k) = fft(abs(real(ifft(output(:,k)))));
    
end

% Half-Hanning FFT * Signal FFT. And then back to time domain
i = [];
for i = 1:numBands
    output(:,i) = real(ifft(sigFreq(:,i).*fft(hann)));
end
